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We present two dimensional numerical simulations of a natural convection problem in an 
unbounded domain. A thermal stratification is applied in the vertical direction and the flow 
circulation is induced by a heat island located on the ground. For this problem, thermal 
perturbations are convected in the horizontal direction far from the heated element so that 
very elongated computational domains have to be used in order to compute accurate numerical 
solutions. To avoid this difficulty thermal sponge layers are added at the vertical boundaries. 
With this approach, stationary solutions at Ra < fO 5 are investigated. Boussinesq equations 
are discretized with a second-order finite volume scheme on a staggered grid combined with a 
second-order projection method for the time integration. 
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In this paper we present numerical simulations of a particular type of thermal fluid flows. Namely, 
we are concerned with the so-called heat island flows 1: i.e. fluid flows where natural convection 
is generated by a local variation of temperature thus inducing buoyancy effect. This phenomenon 
appears in the presence of heat stratification that stabilizes the fluid flow. The present model is 
generally used to study environment problems such as urban heat island [2 [3j . Heat island fluid 
flows occur in open configurations, which require, from the mathematical viewpoint, their study 
in unbounded domains. In practice, numerical simulations are carried out in large but bounded 
computational domains for which appropriate design of boundary conditions must be investigated. 
This issue is a central one in our contribution. 

For this problem, the heat island perturbation generates an ascending flow circulation which 
develops mainly in an area surrounding the heated element. The vertical stratification limits this 
effect by pushing the flow down to the ground. As a consequence of these opposite forces, thermal 
perturbations are propagated in the horizontal direction at long distance far from the heat source. 
Therefore, very elongated domains have to be used in order to accurately compute the temperature 
deviation from the stratified profile. 
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Despite the increase of computational resources, the direct numerical simulation of solutions to 
natural convection induced by a local heat source in large domains remains a real challenge (see [4] 
and the references therein). In most cases, the far- field solutions are unknown so that the use of a 
limited computational area surrounding the heating element requires an appropriate treatment of 
the boundary conditions. In [4], outer artificial conditions are applied at the domain boundaries. 
For heat island flows, only the length of the computational domain has to be limited. Flow 
elevation being strongly reduced by the vertical stratification, the domain length necessary to obtain 
converged solutions is not large. The approach used in this paper consists in applying a thermal 
sponge layer in the vicinity of the vertical boundaries: The temperature equation is modified 
so that the convective terms are smoothly damped in an area closed to the outflow boundaries. 
Boundary conditions with a sponge are classically employed in computational electromagnetics 
[5j [6] and in simulations of compressible turbulent flows [3 [8] . We show that this technique is well 
suited for the numerical simulations of heat island type flows at moderate Rayleigh numbers, that 
is Ra < 10 5 . This issue constitutes the main contribution of this paper. 

The outline of the paper is as follows. In the next section, we describe the set of equations that 
govern the fluid flow in a heat island as well as the domain geometry and boundary conditions. 
In particular, we write the equations in a nondimensional form that involves two parameters, the 
Rayleigh number and a thermal stratification coefficient. We show that a specific choice of relevant 
parameters is to be made here for the heat island flow. Section 3 describes approximations in closed 
domains. Two approaches are used: The former relies on the use of very large computational 
domains without any particular treatment at the domain exits while the latter introduces thermal 
sponge layers acting in the vicinity of the vertical boundaries. Section 4 presents the space and 
time discretization schemes. Preliminary numerical results were obtained by Touzani in [9] using a 
finite element method coupled with a penalty method to impose the incompressibility constraint. 
Here, we use a second-order finite volume scheme on a staggered grid for space discretization and 
a second-order projection method for the time integration of the resulting system of differential 
equations. Section 5 gives numerical results. Numerical simulations in a square differentially 
heated cavity are first performed to check the code accuracy. Stationary solutions to the heat 
island problem at Ra < 10 5 are obtained in very elongated computational domains and are used 
as references to validate the sponge technique applied to the heat equation. Accurate stationary 
solutions are then computed with this approach: Characteristic values are listed for reference. 
Finally, conclusions are drawn and perspectives for future works on this problem are discussed. 



2 Description of the problem 

2.1 The physical problem 

We consider a fluid that fills the half plane {x* = (x* , y*) € K 2 ; y* > 0}. Here and in the sequel 
we shall append a superscript * to all physical dependent and independent variables, the notation 
without * being reserved to nondimensional variables. The fluid is initially at rest and is ther- 
mally stratified in the vertical direction, namely the velocity field u* = (u*,v*) and the potential 
temperature T* satisfy, at time t* = 0, 



where To > is the potential temperature at the ground and a s > is the thermal stratification 
coefficient. 

In order to generate a flow, a local temperature perturbation of intensity T\ > is applied on 



u* = 



T + a s y*, 



(1) 
(2) 
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a source line Q* = (—6/2, 8/2), 8 > 0, located on the ground (see Figure 1), that is we impose 

T^,n = T + ^( 1 -tanh(^l±^)), (3) 

for all i*el 2 n {y* — 0} and for all time t* > 0. Note that this thermal perturbation is constant 
in time and is a regularized version of a heat island type perturbation (see [3] for instance) for 
which a constant and uniform temperature would be applied on the heated element Q*. The 
parameter £ > in ((3]) is used to set the sharpness of the temperature gradient dT*/dx* near the 
plate boundaries \x* \ = 8/2 at the ground level y* — 0. In this study, the value ( = 2.5 x 10~ 2 is 
used. 

Due to the perturbation ([3]), a thermal plume develops above the heated plate Q* . Natural 
convection induces an ascending flow circulation while the gravity force and the vertical stratifica- 
tion limit the development of flow structures in the vertical direction. As a consequence of these 
opposite forces, thermal perturbations are propagated in the horizontal direction at long distance, 
far from the heated element. We may then expect that the solutions will decay rapidly in the 
vertical direction and very slowly in the horizontal one. The main difficulty for this problem, as 
long as numerical simulations are concerned, resides in a suitable design of boundary conditions in 
order to properly reproduce the behavior of the far-field solutions. Indeed, errors in the numerical 
approximation at long distance from the heated plate may deteriorate the accuracy of simulations 
in the region surrounding the heated element. 



2.2 The governing equations 

We consider the set of equations describing a two-dimensional thermal flow assuming Boussinesq 
approximation. Let us recall that this one stipulates that for small temperature differences, the 
density variations are more significant in the gravity acceleration term than in others. 

Velocity u* , pressure p* , density p* and potential temperature T* satisfy the set of equations: 

811* 1 n* 

— - v A*u* + V* ■ (u* ® «•) + -W = - f -ge 2 , (A) 
dt* po Po 

V*-tt*=0, (5) 

^ - k A*T* + V* ■ (u*T*) = 0, (6) 

u*{x*,0) = 0, T*{x*,0)=T + a s y*, (7) 



where the physical constants are the kinematic viscosity v, the thermal conductivity k and the 
modulus of gravity acceleration g. The unit vector in the vertical direction is denoted by e2, 
namely e 2 = (0, 1). 

For these type of flows, we have adopted the approximation that compressibility is expressed 
by a dilatation equation. Therefore, density is related to the temperature variations with respect 
to the reference state (To,po), by the following equation 

p*=po(l-0(T*-T o )), (8) 

where (3 is the thermal expansion coefficient. Equations d4j)~(|8j) are to be considered in the infinite 
domain R x M + where the line y* — contains the heated plate Q* = (—8/2, 8/2). Such a problem 
with appropriate behavior at the infinity is well suited for generating a heat island flow in the 
vicinity of the heated plate if a thermal stratification is given. 

The boundary condition has as effect to generate, in the neighborhood of the plate Q* a 
temperature plume with a shape and an intensity depending on the system parameters. 
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2.3 Nondimensional form of the equations 

In order to normalize the problem, we introduce as usual reference values for temperature T r = T%, 
length L r = 5 and velocity U r = \fgj3L r T r . Reference values for time and pressure can then be 
deduced by t r = L r /U r and p r = po respectively. We define nondimensional space and time 
variables by x — x*/L r and t = t*/t r . In terms of these nondimensional variables, the heated 
plate reads Q = (— 1/2, 1/2). The nondimensional variables 

u* T* - T p* + po g y* 
u=—,T=— and p= (9) 

satisfy inlx R + and for t > the following system of equations: 

3xl i Pt 

— - \ — Au + V- (u®u) + Vp = Te 2 , (10) 
ot V Ra 

V-u = 0, (11) 
dT 1 

:AT + V-(«T)=0, (12) 



dt Vi?a Pr 

u(x,0) = 0, T(x,0) = ay, (13) 
where a — a s L r /T r . The Prandtl and Rayleigh numbers are respectively defined by 



Pr = — and Ra = 



K UK 

The nondimensional form of the boundary condition ([3]) reads 



T(x,t) = Vl-tanh(^+i^ far x = (x,0). (14) 

To simplify, we set the Prandtl number, that characterizes the fluid to its value for the air, Pr = 
0.71. Therefore, the system of equations (jTUJ) — (|T31) depends on two parameters Ra and a. 

The aim of this work is to study through accurate numerical simulations the behavior of the sta- 
tionary solutions and their dependency on the Rayleigh number, at a fixed stratification coefficient 
a. 



3 Approximation in closed domain 
3.1 Large elongated domain 

Clearly, Problem (fTU]) - P^)l is difficult to handle numerically in an unbounded domain. We choose to 
approximate K x R + by a rectangle fl = (-L/2, L/2) x (0, H ) with L large enough (see Figure^ and 
we denote by T its boundary. A simple and naive approach would consist in imposing conditions on 
r consistent with the initial condition (fT3"|) and the boundary condition IT4]) . Therefore, equations 
(fTU)) - (fm]) are supplemented with the boundary conditions: 

u(x,t) = o, ser, (15) 

T(x t t) = Vl-tanh(^^)V xeTo = {xeT; y = 0} 7 (16) 
T(x,t)=ay, ^r\r„, (17) 
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for all time t > 0. 

Such conditions can be enforced only if the domain lengths L and H are large enough so that 
their change has a negligible effect on the solution. As it was previously mentioned, the gravity 
force and the stratified profile (fT"7f tend to push the flow down to the ground, limiting the vertical 
convection. Therefore, the flow variables have a rapid decay with respect to the elevation y, so that 
the domain height H does not need to be too large. The heat island perturbation generates an 
ascending flow circulation which is essentially local. On the other hand, thermal perturbations are 
convected in the horizontal direction at long distance far from the heated source line. Therefore, 
very elongated computational domains have to be considered, that is L ^> H and L ^> 1, in order 
to produce accurate solutions. 



3.2 A truncated temperature equation 

We shall see in the development of this study that the above approach is very consuming in terms of 
computational resources inducing strong limitations on the allowable grid resolution. Indeed, if the 
domain length L is not large enough, artificial boundary layers develop at the domain boundaries 
|x| = L/2 inducing an overestimation of the flow variables. These numerical errors deteriorate 
the accuracy of the solutions in the central area where most of the flow dynamics take place. An 
appropriate design of the behavior of the solution close to these boundaries is necessary in order 
to relax the condition L > 1 on the domain length. To do so, we propose to limit the horizontal 
propagation of the thermal perturbation by damping the convective terms in the temperature 
equation in sponge layers close to the domain exits \x\ = L/2. The nonlinear convection term in 
(TT2)) is multiplied by a filter function 

^ 7 (a0 = e-rvsr) (18) 
where er G (0, 1) and p > 1. This yields the modified heat equation 

'' T ' AT + ^(x)V-(mT) = 0. (19) 



dt ^/Ra P 



r 



For the sake of simplicity, we use the same notation T for the truncated and standard temperature 
respectively solution of (TT9)) and (| 1 2|) . The former corresponds to the choice of 7 = 1 and the 
latter to 7 = 0. 

The filter function ip^ rapidly decays when |x| ~ L/2 whereas ?/> 7 — > 1 in the center of the 
computational domain thus reducing to a classical convective term. The effect of is to introduce 
a sponge layer, close to the vertical boundaries, where the convection of temperature is smoothly 
damped through the outflow. As we will see in Section 5, such a treatment allows to significantly 
reduce the size of the computational domain Q required to reach a given accuracy. 

This approach, while different in its implementation, is similar to the perfectly matched layer 
(PML) method used in computational electromagnetics and introduced by Berenger [5] . Boundary 
conditions with a sponge are also commonly used for the numerical simulations of compressible 
turbulent flows as jet flows for instance [8]. For these problems, the solution is driven to a specified 
outflow state by adding in the sponge layer a cooling term to the right-hand side of the equations. 
We found that our method, for the heat island problem ^0|) - (Tl4|) . is less sensitive to the parameters 
involved in the definition of the sponge functions. 



3.3 A model formulated in terms of temperature fluctuations 

By noting that a stratified profile of the temperature can be expressed in the momentum equation 
as a gradient term, we decompose the potential temperature T into T — ay + 9 which introduces 
the temperature fluctuation 9. By reporting this decomposition into (fTU|) . (fTTj) . (fT9")) and (fT5|) . 
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recalling that u = (u,v) and, introducing the new pressure variable P = p — ay 2 /2, we finally 
obtain the system of equations: 

du I Pt 

— -W — Ati + V -(u®«) + VP = e 2 , (20) 
ot V Ra 

V-it = 0, (21) 

7 --^==A0 + ^ ( (x){V-(ue)+av)=O, (22) 
C£ y/RaPr 

u(x,0)=0, 0(x,O) = O, (23) 
which is supplemented with: 

u(x,t) = 0, ier, (24) 

1 A , /2|x| + 1> 

6(x,t) = o, xer\r„. (26) 



0( !E ,t) = ih-tanh(^^)j, a er„ = {^r ; y = o}, (25) 



In the next section, the numerical approximation of (|20jl - l|25)l is addressed. 

4 Numerical approximation 

The numerical discretization of (|20[) - (|26p is achieved by using a second-order projection scheme in 
time coupled with a second-order finite volume approximation in space. The unknowns are placed 
on a staggered mesh as for the classical MAC scheme [ID] . 

4.1 Time discretization 

The natural convection problem (f2D]) ~ (f!?D|) is solved in two steps decoupling the computation of 
the temperature fluctuation and of the velocity-pressure unknowns. A second-order projection 
scheme [TTJ [HI [T3] is first applied to solve the momentum equations ([2"D)l and to enforce the 
incompressibility constraint (pH]) . 

Let St > stand for the time step and t — kSt discrete time values. Let us consider that 
(vP , Pi, J ) are known for j < k. The computation of (u k+1 , P fe+1 ) consists in: 

- Computing a predictor u k+1 by solving: 

~fc+i t rp, ~k+i . h 1 
u — u K Pr . (u +u K \ 1 



LL A "7" + VP k = ± (3 6 k - e*- 1 ) e 2 



St V Ra V 2 / 2 

- - V • {u k <g> tt fc ) + i V ■ (tt^ 1 (gj-i/- 1 ), (27) 

w fc+1 = on T. (28) 
Projecting to obtain a divergence free velocity u k+1 : 

u k+l _ 1 

? + 77 V(P fe+1 - P k ) = 0, (29) 

St 2 

V • u k+1 = 0, u k+1 ■ n = on L. (30) 
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Finally, the temperature variation 8 k+1 is computed by solving: 

si y/mrp? V 2 J 

= -^(*) (| ( V ' («***) + ~ \ (V • (u*- 1 ^ -1 ) iaw*- 1 )) , (31) 

fc+1 = i(l-tanh(^±±)) on To, fe+1 = on T \ Tq. (32) 

Hence, viscous and diffusion terms are discrctizcd with a Crank-Nicholson scheme while non- 
linear convective terms are integrated by an Adams-Bashforth scheme. The scheme (|27 |) -([32 |) is 
globally second-order accurate. This method is well-suited for the Navier-Stokes equations and is 
frequently used (see for instance [14] and [15]). 

Finally, note that the temperature is computed once the projected (divergence free) velocity is 
obtained. A different approach is used in [TBI H3. HB] : the temperature is first computed, then the 
velocity field is obtained by using the temperature at the new time level. As we will see in Section 
5.1, both time stepping schemes achieve a second-order accuracy and provide similar results. Also, 
in [TB] and [T7], a BDF projection scheme is implemented to discretize the advection-diffusion 
terms. The BDF scheme may be more efficient for the computation of nonstationary solutions as 
it provides a more accurate approximation of the pressure. In the case of stationary solutions, which 
are our main concern in this paper, the more classical approach used here provides satisfactory 
results. 



4.2 Space discretization 

4.2.1 Mesh and unknown locations 

Due to the combined effects of the gravity force and the vertical stratification, flow variables decay 
rapidly with respect to the vertical elevation. Therefore, the domain f2 is discretized by using a 
uniform subdivision in the y-direction. A non-uniform grid is required in the x-direction as L 3> H . 
Let N and M denote two integers and let 

Xi = -if{il) fori = 0,...,N, l=jj, 

V) = jh for j = 0, . . . , M, h = —. 

The function ip : (0, L) — ► (— 1, 1), describing mesh density, is defined by 

2x - L + 7i tanh(7 2 x) - 7i tanh(7 2 (L - a;)) 

V( x > = 1~\ 1 — \Ti — T\ ' ( 33 ) 

L + 7i tann(72-L) 

The function tp enables defining a grid with steps li = x% — a?i-i and distortion ratios r« = £j/^,-_i 
that increase in function of the distance from the center of the heated element. The parameters 
7i and 72 are chosen so that the lengths £i are of order h in the neighborhood of the heated plate 
Q = (-1/2, 1/2). 

We introduce points x i+1 / 2 ■= x%+ ^ +1 for i = 0, . . . , N — 1, and yj+1/2 '■= yj+ | J+1 for j — 
0, . . . , M — 1. All terms in equations (|2T|) . (|2"!)j) . (p3|) and (J2J) are discretized in space by using 
second-order centered finite volume schemes. The discrete unknowns are given on a staggered grid 
(see [10]): discrete pressure values are located at the center of mesh cells K i _i j_i = (xi-i,Xi) x 
(yj-l, Uj), vertical velocity and temperature values are located at the center of mesh cells K { _i j = 
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(xi—i, Xi) x (yj_i, and those of the horizontal velocity are located at the center of mesh 

cells K i j_x = {x i _i 1 x i+ i) x (yj_ 1 ,yj) : as it is shown on Figured We define the vector u k+1 g 
R (JV-i)M f com p 0n ents u^ 1 and similarly, v k+1 g R ^(m-i) jP *+i g r (jv-i)(m-i) and ^fe+i g 
r at(m-i) f components v k +\ p^ 1 and 6^ +1 respectively. 

In [17], the discrete temperature is located at the pressure nodes ,yj_i). However, this 

choice implies interpolations in order to compute the contribution of temperature in the vertical 
velocity momentum equation and vice versa. A more convenient choice for natural convection 
problems is to place temperature at the same nodes as the vertical velocity. 



4.2.2 Discrete systems 



The discretization of (|2T|) is achieved by integration of the equation of the horizontal (resp. vertical) 
velocity component over the volume cells K i j_x (resp. -fQ_i j)- Gradient and Laplace operators 
are classically approximated by centered seconci-order finite volume schemes. Approximation of 
the nonlinear terms requires second-order interpolations of velocity components at the interfaces, 
for instance we use 

CV3 2/ W . f U ij + U i-l,j\ 

u\ Xi _i,y) dy^h^—^ J -j 



Vj-i 



and 



/ 1 

(uv)(x,yj)dx w 



u 



i+hj 



i+l v ij + ^i v i+l,j) 
{ii+l+h) 



Similar interpolation rules are also applied to discretize the equation satisfied by the vertical 
velocity v. This leads to the system of equations 



-fc+i 



6t 
~2 



A^ 1 = -StG!P k 
St 



u k - 



St 



A lU k 



+ 



St 




Pr . ~ k +i 
Ra 



-StG 2 P k + v k 
5 -l(3N 2 (u k ,v k ) 



Ni(u k -\v k - 1 )) 



St 
~2 




A,v k 



N 2 (u k -\v k - 1 )) 



(34) 



(35) 



where the matrices Aj are discrete approximations of the operator — A with appropriate treat- 
ment of the boundary conditions for the velocity components, Gi are the ones of the gradient 
components and Ni are the ones of the nonlinear terms. 

With the use of the staggered implementation of the discrete values on the mesh, several 
possibilities are offered for the treatment of boundary conditions. Concerning the vertical velocity 
component, we choose to impose the boundary conditions on vertical boundaries at grid points 
{{xo,yj),j — 1) • • • , M — 1} and {{xn, yj),j — 1, ■ . ■ , M — 1}. This yields a modified formula for 
the discretization of at the first point away from the vertical boundary, that is 



— (x,y)dx dy w h 



V 3 /2,j ~ Vl/2,j Vl/2,3 ~ U 0J 



x 3/2 ~ x l/2 



%l/2 ~ Xq 



A similar formula applies at the last inner point in the horizontal direction, that is x 



N-l/2- 



On 



horizontal boundaries, boundary conditions for v are imposed at points {(x^_ i , yo), i — 1, ■ ■ ■ , N} 
and {(cCj_ i , yu), i = N}. In the vertical direction, the use of second-order centered formula 
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and of uniform mesh points allows us to apply a discrete Fourier transform [TO] . We thus obtain 
a set of independent and symmetric tridiagonal systems which can be efficiently solved with the 
LDL T algorithm. 

Boundary conditions for u on the vertical boundaries are also imposed at mesh points, that 
is {(xQ,i/j_i) 1 j = 1, ...,M} and {(ccjv, Dj-i )> j — 1, ...,M}. However, at the top and bottom 

horizontal boundaries, values at ghost points J/_i/2 = —\ and Um+i/2 = H + 4 are used to impose 
boundary conditions with a second-order extrapolation formula: We introduce ghost velocity values 



H-i 



= —Uii and UiM+i = ~Ui,M for i = l,...,N. 



The discretization of d 2 /dy 2 on the sequence of mesh points 2/1/2, ■ • ■ , Vm—i/2 with a second-order 
centered finite volume scheme also yields a discrete operator which can be easily diagonalized by 
applying a discrete Fourier transform [19j . 

The discrete version of (|29[) is obtained similarly: 

u k+1 =u k+1 --G 1( f> k+1 , 

2 (36) 

where <f> k+1 = P k+1 — P k . Note that due to the staggered locations of the unknowns, no boundary 
conditions for the pressure are required in the correction step (|36p . Therefore, discrete pressure is 
defined only at interior points. 

The discretization of the incompressibility constraint is achieved by integrating (|3T))) over the 
pressure cell K i _ij_i, leading to 

D lU k+1 + D 2 v k+1 = 0, (37) 

where D\ and D2 are approximations of d/dx and d/dy. Combining (I36[) and (I37| . we deduce the 
linear system satisfied by 0, namely 

(A Gi + D 2 G a ) = + £> 2 « fe+1 ) . (38) 

Once (|3"5)) is solved, the velocity is updated with (|3"r?|) . The linear system defined by (|3"5|) can be 
solved efficiently by applying the same discrete transform used for the vertical velocity component. 
The temperature equation (|3Tj) is integrated over the volume cells K i _i j ; leading to 

Qk+ i St M e u+i = 9 k _ St ^ Qk 

Pr 2y/Ra Pr 

(3N 3 (u k ,v k ,0 k )-N 3 (u k - 1 ,v k -\0 k - 1 j) (39) 

aSt 
~2 



(3^-^), 



where ^ = {^'} e I lV ' M " 1 ', ^ = ip^x^i). Boundary conditions for temperature are treated 
as for the vertical velocity component. 



5 Numerical results 

The main purpose of this paper is to produce reference stationary solutions for heat island flows 
at Rayleigh numbers Ra < 10 5 . In our study the stratification coefficient a is fixed as a = 1. 
Dependency of solutions upon this parameter will be addressed in further works. 
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First, the accuracy of our code is evaluated by computing stationary solutions in a square 
differentially heated cavity. For this test case, benchmark solutions available in [16] are used for 
comparison. Then, stationary solutions for flows in a heat island are studied. The methodology 
used to produce accurate results is detailed and solutions are described and analyzed. 

5.1 Validation of the code: the square differentially heated cavity test 



In order to assess the validity of our code and to check the accuracy of the numerical scheme (|34|) - 
()39l) we have performed numerical simulations of stationary solutions to the square differentially 
heated cavity for values of the Rayleigh number Ra — 10 6 , 10 7 and 10 8 . 

Le Quere [16] produced accurate benchmark solutions for this problem. In [16], Chebyshev 
polynomials were used for the spatial approximation and an influence matrix technique was applied 
in order to enforce the divergence free condition. Note that this problem, described in [20l [Pol fTB] . 
differs from the heat island problem by the computational domain and the boundary conditions. 
However, the discrete system (|34l) - (|3l?|) . with 7 = and a = 0, applies as well to this test case. 

Stationary solutions were obtained on uniform grids, in both horizontal and vertical directions, 
with mesh sizes decreasing from 1/32 through 1/1024. The choice of uniform grids is not optimal 
for this problem as boundary layers develop along the vertical heated walls: A large number of 
points is thus required in order to obtain accurate results. Such a choice is however convenient 
and allows us to easily check the code accuracy. 

The characteristic values suggested by De Vahl Davis in [5D] were computed and compared with 
those of the benchmark solutions [16]. All these values are recovered and a second-order spatial 
convergence is obtained (see Figured]). 

5.2 Stationary solutions of flows in a heat island 

Due to the presence of the vertical stratification, the thermal perturbations are convected in the 
horizontal direction far from the heated source line. As a consequence, very long domains have 
to be used in order to produce accurate results. Numerical simulations in small computational 
domains are contaminated by artificial boundary layers which develop at the outflow boundaries 
\x\ = L/2. If the domain length L is not large enough, the temperature cannot smoothly relax 
towards the vertical stratified profile imposed on the boundaries. 

For fixed Rayleigh numbers and mesh sizes, numerical simulations in domains with increasing 
lengths are performed with the standard heat equation, that is pip with 7 = 0. The effects of L 
and H on the accuracy of the results are investigated. This approach while time consuming allows 
to produce reference solutions. Numerical simulations with the truncated temperature equation 
(7 = 1) are then performed for comparison. This study demonstrates the efficiency of the thermal 
sponge layers. Finally, stationary solutions at Rayleigh numbers Ra < 10 5 are computed on meshes 
with a finer resolution. Characteristic values are reported and various profiles are reproduced and 
analyzed. 

The stationary state of the numerical simulations was assumed to be reached when time vari- 
ations of flow variables arc controlled as it follows 



where Tol £ (10 10 , 10 8 ) is a given parameter. 

5.2.1 Numerical simulations in large elongated domains 

For fixed vertical resolutions, h = 1/16, 1/32 and 1/64, and Rayleigh numbers Ra = 10 3 , 10 4 
and 10 5 , numerical simulations were performed for increasing values of the sizes L and H of the 



case 
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computational domain. The discrete systems ([3~4")) - (f3T?| with the parameter 7 set to zero was used 
which corresponds to the classical heat equation. The temperature fluctuation 9 is the flow variable 
for which convergence with respect to the domain size is the slowest. Therefore, we choose to use as 
reference value to compare simulations in different domains the minimum value reached by 9 inside 
the computational domain. For stationary solutions, the minimum is reached above the center of 
the heated element, namely on the axis x — and for y > 0. Recalling that fi = (-§, §) x (0, H), 

we denote = mint Xt y) e n 6(x, y) and we define the reference values (L Te {, H IC [) and O^in by 

\n{L,H) _ grei I 

HLM) ■= min | flrc f 1 mi < 0- 1 h fOT L > L ^ H > H ^ (4°) 

Fminl 

so that 9™f n is considered as a converged value for fl^jj^ . Values (L re {, H rc {) found for the consid- 
ered Rayleigh numbers and mesh sizes are reported in Table I. We observe that the domain length 
L re f is not sensitive to the Rayleigh number while the domain height H Ie f, for a fixed resolution, 
decreases when Ra increases. The strength of the stratification grows with the Rayleigh number: 
The flow is pushed down to the ground. 

Once these reference domain sizes have been found by repeated simulations in large domains, 
values of (L, H) ensuring a h 2 approximation of 9^'^ were estimated. We introduce (L c , H c ) so 
that ei y L a ,H a ) < h 2 an d 

e( L ,H) > h 2 for L<L C ,H< H c . (41) 

Therefore, for fixed Rayleigh number and mesh size, (L c , H c ) are the minimum values of the domain 
sizes required to compute a numerical solution accurate up to scheme accuracy. Estimates of these 
minimal values are reported in Table I. The tests (j4"0|) and (|4"Tj) impose strong restrictions on the 
admissible domain sizes. Indeed, errors are most often larger than h 2 even for second-order schemes 
(see Figure S] for example). The values (L C ,H C ) estimated with (|4"Tj) are probably too restrictive. 
However, their use ensures accurate results. 

Figures \5\ represents the convergence history of £(l,h) with respect to the domain length L for 
H = 4, 6 and 8. The Rayleigh number is Ra = 10 4 and the vertical resolution is h — 1/32. The 
convergence rate of towards the reference value 0JJjf n behaves like l/L for small values of 

L and like L~ 3 / 2 for large values of L. The scheme accuracy is reached for H — 6 and L of the 
order of 2 000 (see also Table I). The same behavior with respect to L was found for other vertical 
resolutions and Rayleigh numbers. 

On Figure [H the convergence history of C(l,h) f° r H — H c is represented for Ra — 10 3 , 10 4 
and 10 5 . The convergence rate of (-(l.h) is found to be independent of the Rayleigh number: All 
curves have the same slope in logarithmic scales. Also, the value h 2 = 1/32 2 corresponding to the 
scheme accuracy is reached by £(l,h c ) f° r values of L decreasing when Ra is increased. 

For the vertical resolution h = 1/64, the computation of the reference solution in the domain 
f^ref = (0, £ re f) x (0, -ff r ef) with L re f = 10 000 was achieved with 28 000 points in the horizontal 
direction. This guarantees that the mesh satisfies li w h in the neighborhood of the heated element. 
With the use of the mapping function (f3"3")l , we have in that case: max^/ min£i = 74. Therefore 
35.8 (resp. 17.9) millions of points were used to compute the reference solution at Ra = 10 3 (resp. 
10 5 ) and 35 000 (resp. 200 000) time iterations were required to reach the stationary solution. This 
represents 600 (resp. 1 500) monoprocessor computing hours on an IBM Power4 computer. 

Due to this need for a large amount of computational resources, such a study is not feasible on 
meshes with smaller grid sizes in the vertical direction. By extrapolating the values L c and H c found 
for h — 1/16, 1/32 and 1/64 (see Table I), we can roughly estimate these minimal values for a finer 
mesh, namely h = 1/128. Wc obtain {L C ,H C ) = (12 000,18) at Ra = 10 3 , (L C ,H C ) = (10 000, 12) at 
Ra = 10 4 and (L C ,H C ) = (8 000, 10) at Ra = 10 5 . Therefore, a resolution of at least 115 millions 
of points would be required in order to compute a reference solution on a grid with a vertical 
resolution h = 1/128. As it is shown in the next sections, the use of the truncated heat equation 
(7 = 1) allows us to significantly relax these constraints on the computational parameters. 
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5.2.2 Efficiency of the truncated temperature equation 

The effect of the filter function (|18|) is to smoothly damp the convective terms in the heat equation 
in the vicinity of the domain boundaries |x| = L/2. By introducing such thermal sponge layers, 
we aim at improving the accuracy of the numerical simulations when the computational domains 
are not long enough to ensure a 0(h ) approximation. 

Numerical simulations of stationary solutions at Ra = 10 5 have been performed for H fixed to 
H c , for the vertical resolutions listed in Table I and for increasing values of L. The errors e(L, H c ) 
produced by the standard and the modified heat equation are used to compare the efficiency of 
both models. The values a — 0.85 and p = 8.0 have been retained for the filter function (fT8|) . 
These parameters were found to be efficient for the numerical simulations of stationary solutions. 
A parametric study of the truncated temperature equation is beyond the scope of this paper. 
However, this question is important and will be addressed in future works on this problem. 

On Figure [3 the errors e(L,H c ) obtained for h = 1/32 with the standard and the truncated 
heat equations are plotted. We note that: 

— the truncated heat equation produces errors about 10 times smaller than the standard equation 

even for small values of L; 

— both curves have the same decay rate and converge to the same asymptotic value; 

— the truncated heat equation produces values of L c which are approximately 3 times smaller 

than the values listed in Table I and corresponding to the standard heat equation. We recall 
that L c , defined by (|4ip. is the minimum value of the domain length required to ensure a h 2 
approximation of the temperature fluctuation. 

The values L c = 120, 400 and 900 are obtained with the truncated equation for the respective 
resolutions h = 1/16, 1/32 and 1/64 while L c — 480, 1280 and 3 200 were necessary with the 
classical heat equation (see Table I). 

Also, by examining the time history of the discrete time variation - ^ — — , it appears that 

the convergence to the stationary solution is achieved in less time iterations with the truncated 
equation than with the classical one. For example, in £1 = (—240,240) x (0,4) and for h = 1/32 
the stationary solution at Ra = 10 5 is reached after 32 000 time iterations with the truncated heat 
equation while 48 000 time iterations are needed with the classical one (see Figure [5]) . 

Therefore, for a given accuracy, stationary solutions can be computed with the truncated 
heat equation in significantly smaller computational domains than with the classical temperature 
equation and in less time iterations. This results in a use of less computational resources. As 
a consequence, this approach allows us to compute stationary solutions on meshes with a finer 
vertical resolution. 

5.2.3 Accurate stationary solutions 

Direct numerical simulations at Ra = 10 3 , 10 4 and 10 5 and with a vertical resolution h — 1/128 
have been performed with the truncated temperature equation. The computational parameters are 
listed in Table II. The estimates derived in Section 5.2.1 are used to determine the computational 
domains. Also, in agreement with the previous section, domain lengths about 3 times smaller than 
the estimated values are retained. 

In order to characterize the stationary solutions, the maximum values of the velocity compo- 
nents (it, v), the temperature variation 9, the vorticity uj = dv/dx — du/dy and the streamfunction 
ip are reported in Table III. The locations in the computational domains where these extrema are 
reached are also collected: When only one of the coordinates is listed, the other one is equal to 0. 

Also reported in Table III is the Nussclt number which is defined by 
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As expected, this value, which measures the intensity of the heat transfer, increases with the 
Raylcigh number. 

The velocity components and the temperature fluctuation decay rapidly with respect to the 
elevation y as it is shown on Figures [9] and [11] representing the vertical profiles at the center of 
the heated element x — 0. The gravitational force and the vertical stratification limit the vertical 
propagation of perturbations. In the horizontal direction (see Figure [TUb), the vertical velocity 
v vanishes rapidly for x outside of the heated region, that is for \x\ > 0.5. Therefore, vertical 
convection is essentially localized above the heated element: This behavior is independent of the 
Rayleigh number. However, its intensity increases with Ra. Indeed, the maximum value reached by 
the vertical velocity component increases with Ra (see Figure [TUb and Table III). The temperature 
fluctuation and the horizontal velocity have a similar behavior in the horizontal direction for \x\ > 2 
(see Figures [TUb and [P2f : They decay slowly to a small but nonzero value which is growing with 
Ra. Therefore, the convection outside the heated source line is mainly horizontal. This illustrates 
the difficulty to approximate such flows in limited computational domains. 

The profiles of the temperature fluctuation, plotted on Figures and [POk . show that the as- 
cending propagation of the thermal perturbation is reduced when Ra is increased. Simultaneously, 
the profile of the velocity components exhibit largest extrema and steepest gradients for y < 1 
(see Figures [9]d and fTTj) . Also, elevations where the velocity components are maximum decrease 
for growing Ra (see Table III). Hence, when the Rayleigh number is increased, the flow is pushed 
down to the ground. At larger Rayleigh numbers, we expect that the competition between the 
natural convection, inducing an ascending propagation, and the vertical stratification, limiting this 
effect, will induce a loss of symmetry of the solutions leading to unsteady flows. 

To better illustrate the effect of the vertical stratification, isolines of the temperature fluctuation 
9 and the vorticity u> are displayed on Figures [13] and [14] in a region surrounding the heat island 
perturbation, that is for \x\ < 5 and y < 3. The thermal plume in form of a mushroom, typical in 
natural convection problems (see [3] for instance), cannot develop in a stratified medium. Instead, 
the main thermal structure is centered above the heated plate, symmetric with respect to the axis 
x = and very elongated in the horizontal direction. Above, a thermal sink characterized by 
negative temperature variation 9 is observed. The intensity of this structure grows with Ra while 
its vertical position decreases. The vorticity structures (see Figure [T4|) exhibit multi-cell symmetric 
patterns. They become thinner when Ra increases and are clearly pushed down to the axis y = 0. 

5.2.4 Computational efficiency 

A parallel version of the Fortran 90 code based on implicit communications (OpenMP) was used for 
the numerical simulations presented in this paper: An efficiency of approximately 6.8 is found on 8 
processors on a cluster of IBM Power 4 computers. In order to perform the numerical simulations 
presented in Section 5.2.3, 6 000 monoprocessor hours were necessary. The CPU time per iteration 
and per node used by the code is 2.5 x 10 -6 seconds on IBM Power 4 processors. Concerning the 
memory, 20 real unknowns (8 bytes) have to be stored for each node of the mesh. 

6 Concluding remarks and perspectives 

In this paper, steady state solutions of a natural convection problem in an unbounded domain are 
investigated by direct numerical simulations. For this problem, the flow is thermally stratified in 
the vertical direction and perturbed by a local heat island located on the ground. Due to the verti- 
cal stratification, the flow circulation is dominated by horizontal convection, so that perturbations 
are propagated in the horizontal direction far from the heated source. Stationary solutions are first 
investigated by numerical simulations in very elongated domains for moderate vertical resolutions, 
that is h = 1/16,1/32 and 1/64. Repeated computations in increasing domains have been per- 
formed: the minimum length and height necessary to ensure a 0{h?) accuracy have been estimated 
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at Ra — 10 3 , 10 4 and 10 5 . This approach, while time and memory consuming, provided reference 
simulations that have been used to validate and compare results obtained with a truncated heat 
equation. We have shown that the use of a suitable thermal sponge layer placed at the vertical 
outflow allows to noticeably reduce the size of the computational domain. Therefore, numerical 
simulations on finer grids are made accessible. The stationary solutions at the aforementioned Ra 
have been computed on grids with vertical resolution h — 1/128. Characteristic values of these 
steady states have been provided. 

The thermal circulation induced by the heat island consists in symmetric multi-cell pattern 
centered above the heated element. Flow structures are pushed down to the ground when the 
Rayleigh number is increased. Also, their intensity grows with Ra. We therefore may expect 
that stability of steady states will be lost at larger Ra leading to nonstationary solutions. The 
thermal sink found above the heat island should first oscillate with respect to the vertical axis 
x = in a periodic time regime. The numerical study of the development of instabilities and the 
detection of successive transitions from steady state to turbulent flows is our main motivation. 
Contributions to this project will be presented in forthcoming papers. Dependency of solutions 
upon the stratification coefficient is also an open question for this problem. Such study will be 
addressed in future works. 
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Figure 1: Heat island perturbation. 
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Figure 2: Computational domain Q = (— -§,■§) x (0,£f). 
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Figure 3: Cells K i _i j_i (solid), K^i ■ (dashed) and K i -_i (dotted) and their corresponding 
discrete values. 
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Figure 4: Spatial accuracy for the square differentially heated cavity test case at Ra — 10 6 . Solid 
line: slope 2 ; dashed line: horizontal velocity ; dotted line: vertical velocity ; dashed-dotted line: 
Nusselt number. 



Table 1: Estimate of the size of the computational domain for various mesh size h and Rayleigh 
number Ra. 



Ra 


Ny/H 


(L C ,H C ) 


(Lref , -Href) 


10 3 


16 


(960,8) 


(3 200, 12) 




32 


(2 560, 10) 


(6 200, 16) 




64 


(5 200, 14) 


(10 000,20) 


10 4 


16 


(640, 5) 


(3 200,8) 




32 


(1920,6) 


(6 200, 10) 




64 


(4 200,8) 


(10 000,12) 


10 5 


16 


(480, 3) 


(3 200,6) 




32 


(1280,4) 


(6 200,8) 




64 


(3 200,6) 


(10 000,10) 



Table 2: Computational parameters for numerical simulations on meshes with h = 1/128. 



Ra 


10' 5 


10" 


10 s 


(L,H) 


(4400, 18) 


(3 600,12) 


(2 560, 10) 


(N , M) 


(19 000,2 304) 


(16 500,1536) 


(12 000, 1280) 


St 


0.1 


0.05 


0.025 


T stat 


4 893 


6 569 


10 902 
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Figure 5: Convergence of £(l.h) at Ra = 10 4 with respect to the domain length for H = 4 (dotted 
line), H = 6 (dotted-dashcd line) and H = 8 (solid line). The dashed line corresponds to the 
expected accuracy h 2 = 1/32 2 . 
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Figure 6: Convergence of £(l,h c ) with respect to the domain length for Ra — 10 3 (solid line), 
Ra = 10 4 (dotted line) and Ra — 10 5 (dashed-dotted line). The dashed line corresponds to the 
expected accuracy h 2 = 1/32 2 . 
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Figure 7: Convergence of £(l.h=4) at Ra — 10 5 with respect to the domain length obtained with 
the standard heat equation (solid line) and the truncated version (dashed-dotted line) . The dashed 
line corresponds to the expected accuracy h 2 — 1/32 2 . 
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Figure 8: Convergence history of maxi.j j^- 

(dotted line) heat equation at Ra = 10 5 . The vertical resolution is h = 1/32, the time step is 
St = 0.1 and the computational domain is O = (-240, 240) x (0, 4). 
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Table 3: Characteristic values of the stationary solutions. 



Ra 




10 d 




10 4 




10 b 


$min 




-0.024823 




-0.071289 




-0.166316 


V 




1.36191 




0.94640 




0.84320 


Wmax 




0.118887 




0.174844 




0.179054 


{x,y) 


(- 


0.52314,0.17705) 




-0.36781,0.13255) 


(- 


-0.30643, 0.09322) 


"max 




0.125594 




0.228250 




0.322483 


2/ 




0.55124 




0.44684 




0.42552 


^min 




-0.030470 




-0.039291 




-0.079265 


{x,y) 


(" 


■0.85107,0.41913) 


(■ 


-0.60361,0.32122) 


(" 


-0.39114,0.64777) 


^max 




2.06900 




3.951325 




5.921375 


a: 




0.49642 




0.48402 




0.47813 


V'max 




0.042954 




0.048265 




0.040272 


(z,2/) 


(- 


-0.60014, 0.59882) 




-0.38883, 0.45575) 




-0.25610, 0.45867) 


Nu 




0.148605 




0.295132 




0.643594 




Elevation y 



0.35 




Elevation y 



Figure 9: Profiles of the temperature variation 9 (a) and of the vertical velocity v (b) at the center 
of the heated element, i.e. at x = 0, for Ra = 10 3 (solid line), Ra = 10 4 (dashed line) and 
Ra = 10 5 (dashed-dotted line). The vertical resolution h = 1/128 and the truncated heat equation 
were used. 




Figure 10: Profiles of the temperature variation 6 (a) and of the vertical velocity v (b) at the 
elevation y = 0.5, for Ra = 10 3 (solid line), Ra = 10 4 (dashed line) and Ra — 10 5 (dashed-dotted 
line). The vertical resolution h = 1/128 and the truncated heat equation were used. 
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Figure 11: Profile of the horizontal velocity u at x — 0.25 for Ra = 10 3 (solid line), Ra = 10 4 
(dotted line) and Ra = 10 5 (dashed line). The vertical resolution h = 1/128 and the truncated 
heat equation were used. 




Figure 12: Profile of the horizontal velocity u at the elevation y = 0.1 for Ra = 10 3 (solid line), 
Ra = 10 4 (dotted line) and Ra — 10 5 (dashed line). The vertical resolution h = 1/128 and the 
truncated heat equation were used. 
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Ra = 10 3 



Dashed lines: 
-O.OOl. 



x 

-0.024, -0.02, -0.015, -0.01, -0.008, -0.007, -0.006, -0.005, -0.004, -0.003, -0.002, 




Solid lines: 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.078, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 1.0. 
Dashed lines are opposite values 



